Does the Brain’s E/I Balance Really Shape Long-Range Temporal Correlations?

Lydia Sochan; Alexander Mark Weber

¹ School of Biomedical Engineering, The University of British Columbia, Vancouver, BC, Canada
² BC Children’s Hospital Research Institute, The University of British Columbia, Vancouver, BC, Canada
³ Pediatrics, The University of British Columbia, Vancouver, BC, Canada

^✉ Correspondence: Alexander Mark Weber <aweber@bcchr.ca>

In [1]:

# options
knitr::opts_chunk$set(
  # fig.width=8, fig.height=5,
  # out.width="50%",
  # fig.align="center",
  echo=FALSE,
  message=FALSE,
  warning=FALSE
  # cache=TRUE
)
set.seed(1234) # reproducible
options(knitr.kable.NA = '') # how kable handles NA
options(reticulate.repl.quiet = TRUE)

In [2]:

#libraries
library(tidyverse) # ggplot2, dplyr, tidyr, readr, purrr, tibble, stringr, forcats

── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.4     ✔ readr     2.1.5
✔ forcats   1.0.0     ✔ stringr   1.5.1
✔ ggplot2   3.5.1     ✔ tibble    3.2.1
✔ lubridate 1.9.3     ✔ tidyr     1.3.1
✔ purrr     1.0.4     
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

library(english) # numbers to words
library(tools) # title-case
capitalize_first <- function(x) {
  paste0(toupper(substr(x, 1, 1)), substr(x, 2, nchar(x)))
} # capitalize first letter
library(patchwork) # figure manipulation
library(magick) # load images

Linking to ImageMagick 6.9.11.60
Enabled features: fontconfig, freetype, fftw, heic, lcms, pango, webp, x11
Disabled features: cairo, ghostscript, raw, rsvg
Using 32 threads

library(ggpubr)

In [3]:

load("notebooks/results.RData")

In [4]:

format_p <- function(p, digits = 2) {
  if (p < 0.01) {
    "< 0.01"
  } else {
    format(round(p, digits), nsmall = digits)
  }
}

Abstract

1 Introduction

Thirty years ago, functional magnetic resonance imaging (fMRI) profoundly changed the world of MRI by allowing real-time analysis of pressing neuropsychological questions^1–3. While initially used to probe task-based responses, researchers have more recently developed an interest in studying brain function at rest, known as resting-state fMRI (rs-fMRI)⁴, i.e. to understand how brain dynamics at rest are related to neurological functioning as well as individual differences. A critical tool in analyzing these dynamics is the Hurst exponent (H)⁵, a measure of self-similarity derived from the blood-oxygen-dependent (BOLD) signal. H estimates the extent to which the BOLD signal displays long-term memory, where a higher value indicates a self-similar signal with long-term positive autocorrelations^5,6. Another way of understanding H is that a signal with high H is fractal: similar temporal signal fluctuations are observed, no matter the time scale⁵.

H has also emerged as a valuable tool in clinical research, capturing changes in BOLD signal dynamics across various neuropsychiatric conditions. In aging populations for instance, H has been found to be elevated^7,8; this relationship has also been found in mild cognitive impairment and Alzheimer’s disease^9,10. Additionally, changes in H have been observed in conditions such as autism, mood disorders, and brain injury^11–15. These differences suggest H may reflect changes in global and local functioning.

Underlying these observations is the critical brain hypothesis, which posits that the brain operates at a critical point, a state where order and disorder are perfectly balanced to enable optimal information processing^4,6,16–27. When operating at a critical point, the brain is maximally sensitive to external stimuli, and can dynamically transition between ordered and disordered states to promote efficient cognitive processing^4,6,26,27.Recent papers suggest the control parameter underlying the brain’s ability to transition between states is the excitatory-inhibitory (E/I) ratio, the balance of excitatory and inhibitory neural activity, often operationalized as the ratio of the primary excitatory-to-inhibitory neurotransmitters, i.e. glutamate-to-GABA ratio^19,21–25. It is thought that E/I controls criticality by modulating the brain’s signal-to-noise ratio^19,28.

Besides the implications to the critical brain hypothesis, understanding if E/I is related to H may allow for easier estimation of excitatory and inhibitory neurotransmitters, as accurate E/I measurement is technically challenging²⁹. Magnetic resonance spectroscopy (MRS) is the only non-invasive method of measuring the ratio of glutamate (Glu; excitatory) to Gamma-aminobutyric acid (GABA; inhibitory) in vivo³⁰. Unfortunately, it has both limited spatial and temporal resolution^25,29,30. Consequently, if H could serve as a proxy for E/I, it would be much easier to estimate E/I in conditions of interest such as autism, Alzheimer’s, and schizophrenia.

There are a handful of studies suggesting a link between H and E/I, however they are all either computational models or animal studies^19–25. Moreover, their findings are inconsistent, with some reporting positive linear, negative linear, or U-shaped relationships between the two variables (see Table 1). Thus, there is a need for further study, both to clarify the nature of a potential E/I-Hurst relationship, and also to confirm if this relationship indeed exists in the human brain. Therefore, the current study seeks to investigate the potential E/I-Hurst relationship in vivo, within the visual cortex during movie-watching and rest.

Table 1. Summary of methods for existing E/I-Hurst studies

Citation	Study Type	H Data Type	H Calculation Method	E/I Type	E/I-Hurst Relationship
Poil et al. (2012)²⁰	Computational with in-house simulated model	Neuronal avalanche size	Detrendend fluctuation analysis (DFA)	Structural: number of E-to-I neurons	Inverse U
Bruining et al. (2020)²³	Computational with model by Poil et al. (2012); modified in-house	Neuronal oscillation amplitude	DFA	Structural: number of E-to-I synapses	Inverse U
Gao et al. (2017)²⁵	Computational; in vivo in rats and macaques	Local field potential (LFP)	PSD	Estimated from LFP	Positive linear
Lombardi et al. (2017)²¹	Computational with in-house model	Neuronal avalanche size	PSD	Structural: number of E-to-I neurons	Negative linear
Trakoshis et al. (2020)²⁴	Computational with simulated data; in vivo in mice	fMRI BOLD signal	Wavelet-based maximum likelihood method	E-to-I synaptic conductance	Positive linear

2 Methods

2.1 Participants

Twenty-seven healthy adult participants were recruited to the study. One participant was not scanned due to feelings of claustrophobia while in the scanner. After our analysis and performing quality assurance (see below), a further seven participants were removed for having less than ideal MRI data quality, leaving nineteen final participants, between the ages of 21.3 and 53.4 (mean age ± sd: 30.1 ± 8.7 years; 9 males).

2.2 Ethics Statement

Written informed consent was obtained from all participants. Ethics approval was granted by the Clinical Research Ethics Board at the University of British Columbia and BC Children’s & Women’s Hospital (H21-02686).

2.3 Scanning Procedure

After two anatomical sequences were acquired, participants were instructed to visually fixate on a cross-hair for 24 minutes. During this period, an fMRI, MEGA-PRESS, and sLASER sequence were acquired (see Figure 1 A and Section 2.4). Next, participants were instructed to watch a nature documetnary (Our Planet (2019), Episode 3, “Jungles”³¹) for 24 minutes. During this period, another set of fMRI, sLASER, and MEGA-PRESS sequences were acquired. See Figure 1 B for a visual representation of the scanning protocol. Total scan duration was approximately 1 hour. All participants followed the same order of rest than movie, and all participants saw the same movie segment, beginning at the same time during the scan.

In [5]:

fmrimrs <- image_read('../Figures/StudyDesign/studydesign_fmri_mrs.png')
fmrimrs_dim <- image_info(fmrimrs)

fmrimrs <- ggplot() + background_image(fmrimrs) + theme_void() +
  theme(aspect.ratio = fmrimrs_dim$height / fmrimrs_dim$width)
studydesign <- image_read('../Figures/StudyDesign/StudyDesign.png')
studydesign <- ggplot() + background_image(studydesign) + theme_void()
p <- (fmrimrs + plot_spacer()) / studydesign
p + plot_annotation(tag_levels = 'A') & theme(plot.tag = element_text(size = 11))

**Figure** 1. **Overview of the MRI acquisition protocol.** A) fMRI coverage was across the whole-brain (example coverage in red). MRS voxels were placed in the visual cortex (blue). Background image is of a T1w acquisition from a sample subject. B) fMRI, MEGA-PRESS and sLASER were acquired first with participants looking at a white cross for 24 mintues. Next, the same sequences were acquired with the participants viewing a nature documentary for an identical period of time.

2.4 Acquisition Details

Scans were performed at BC Children’s Hospital MRI Research Facility on a 3.0 Tesla GE Discovery MR750 scanner (scanner software version: DV26.0_R03) with a Nova Medical 32 channel head coil. Participants changed into scrubs and were screened by an MRI technologist. Participants were given wax earplugs and a fiducial was placed on their left temple. Participants were provided with an audio headset and blanket once lying down on the scanner bed. Since visual stimuli were to be rear-projected, position and angle of mirror above patient eyes was adjusted for optimal movie viewing.

The following MRI scans were acquired. A 3D T1-weighted sagittal fast spoiled gradient echo (FSPGR) sequence; a 3D T2-weighted sagittal CUBE; a 2D echo-planar imaging (EPI) multi-echo gradient-echo fMRI sequence; a semi-LASER (sLASER) sequence; and a MEGA-PRESS sequence. Details are listed in Table 2.

Table 2. Summary of MRI acquisition details

Sequence	TE (ms)	TR (ms)	Flip Angle	FOV	Slice Thickness (mm)	In-Plane Resolution (mm²)	Other Parameters
3D T1-weighted FSPGR	2.176	7.216	12	256 × 256	0.9	0.9375 × 0.9375	-
3D T2-weighted CUBE	75.242	2,504.0	90	256 x 256	0.9	0.9375 × 0.9375	-
2D EPI Multi-Echo fMRI	12.2, 35.352, 58.504	1,500	52	64 x 64	3.6	3.5938 × 3.5938	Acceleration factor = 6
sLASER	35	2,000	90	-	28.0	28.0 x 28.0	-
MEGA-PRESS	68	1,800	90	-	28.0	28.0 x 28.0	-

For both MRS sequences, the voxel size was set to 2.8 x 2.8 x 2.8 cm³. MRS voxels were rotated and placed in the occipital lobe, aligned along the calcarine fissure. Finally, blip-up and blip-down spin-echo versions of the fMRI sequence were acquired at the end to estimate the B0 non-uniformity map for fMRI phase distortion correction.

2.5 Image Processing

Our full image processing pipeline has been can be accessed from our Github account page: github.com/WeberLab/EI_Hurst_Analysis

Images were downloaded offline from the scanner in raw Digital Imaging and Communications in Medicine (DICOM) format. DICOM files were then converted to Neuroimaging Informatics Technology Initiative (NIfTI) using Chris Rorden’s dcm2niix³² (v1.0.20211006) and then to Brain Imaging Data Structure (BIDS)³³ format using dcm2bids³⁴ (v2.1.6).

2.5.1 Structural Images

The T1w image was corrected for intensity non-uniformity with N4BiasFieldCorrection³⁵ and distributed using ANTs³⁶ (v2.3.335) to be used as a T1w-reference for the rest of the workflow. The T1w-reference was skull-stripped using a Nipype³⁷ implementation of antsBrainExtraction.sh from ANTs; OASIS30ANTs was used as a target template. Fast³⁸ (FSL³⁹ v.6.0.5.1:57b01774, RRID: SCR_002823) was used for brain tissue segmentation into cerebrospinal fluid (CSF), white matter (WM), and gray matter (GM). Brain surfaces were reconstructed with recon-all⁴⁰ (FreeSurfer⁴⁰ 7.3.2, RRID: SCR_001847). The previously-estimated brain mask was refined with Mindboggle⁴¹ (RRID:SCR_002438) to reconcile ANTs-derived and FreeSurfer-derived segmentations of cortical GM. AntsRegistration³⁶ (ANTs 2.3.3) was used to perform volume-based spatial normalization to two standard spaces: MNI152NLin2009cAsym and MNI152NLin6Asym. Normalization used brain-extracted versions of both T1w reference and T1w template.

2.5.2 fMRI

Using fMRIPrep⁴², the shortest echo of the BOLD run was used to generate a reference volume (both skull-stripped and skull-included). Head-motion parameters with respect to the BOLD reference (transformation matrices as well as six corresponding rotation and translation parameters) were estimated before spatiotemporal filtering using mcflirt⁴³ (FSL v6.0.5.1:57b01774). The fieldmap was aligned with rigid registration to the target EPI reference run. Field coefficients were mapped to the reference EPI using the transform. BOLD runs were slice-time corrected to 643 ms (half of slice acquisitio range of 0-1290 ms) using 3dTshift from AFNI⁴⁴ (RRIS: SCR_005927). To estimate T2* map from preprocessed EPI echoes, a voxel-wise fitting was performed by fitting the maximum number of echoes with reliable echoes in a particular voxel to a monoexponential signal decay model with nonlinear regression. Initial values were T2/S0 estimates from a log-linear regression fit. This calculated T2 map was then used to optimally combine preprocessed BOLD across echoes using the method by Posse et al. (1999)⁴⁵. The generated BOLD reference was then co-registered (6 degrees of freedom) to the T1w reference with bbregister (FreeSurfer⁴⁰) using boundary-based registration. First, a reference volume and its skull-stripped equivalent were generated with fMRIPrep. Confounding time series were calculated from preprocessed BOLD: framewise displacement (FD), DVARS, and three region-wise global signals. Tedana⁴⁶ was then used to denoise the data by decomposing the multi-echo BOLD data via principal component analysis (PCA) and independent component analysis (ICA). The resulting components are automatically analyzed to determine whether they are TE-dependent or -independent. TE-dependent components were classified as BOLD, while TE-independent components were classified as non-BOLD and were discarded as part of data cleaning.

Participants were excluded from further ananlysis if their mean FD was > 0.15 mm.

2.6 MRS

sLASER data were processed and fit to a spectrum using Osprey⁴⁷ (v2.4.0). Full width half-maximum (FWHM) of the single-Lorentzian fit of the N-acetylaspartate (NAA) peak were calculated for quality assurance purposes. The MRS voxel was co-registered to T1w reference image and segmented by SPM12⁴⁸ into CSF, GM, and WM. Metabolites were water-scaled as well as tissue- and relaxation-corrected by the Gasparovic et al. (2006) method⁴⁹. Processed (non-fitted) data from Osprey was also fed to LCModel⁵⁰ (v6.1) for fitting. Cramer-Rao lower bounds (CRLBs) were calculated from LCModel fit. Osprey estimates were used for further analysis. Glutamate is challenging to capture due to its signal overlaps with other metabolites⁵¹. In particular, Glu shares a similar chemical structure with glutamine (Gln) which causes the spectral features of Glu to be contaminated with Gln⁵². As a result, we decided to report Glx values, which co-reports Glu and Gln to avoid errors in spectral assignment, especially since it is controversial whether Glu can reliably be separated from Gln at 3T^51,52. We henceforth refer to Glu as either Glx or glutamate.

MEGA-PRESS data were processed and fit with Osprey and LCModel as previously described (see Figure 5). Data were also fit with Gannet (v3.3.0); they were relaxation-, tissue- and alpha-corrected using the Harris et al. (2015) method⁵³. Osprey values were ultimately used for further analysis. Due to the J-editing sequence of MEGA-PRESS, a challenge of GABA quantification is macromolecule quantification⁵³. As a result, we report GABA+, a measure which co-reports GABA with macromolecules. Macromolecules (MM) are expected to account for approximately 45% of the GABA+ signal⁵³. While a macromolecule-suppressed estimate of GABA seems ideal, a recent 25-site and multi-vendor study conducted at 3T found that GABA+ showed much lower coefficient of variation than MM-suppressed GABA, meaning that GABA+ is more consistent across research sites and MRI vendors (i.e., Philips, GE, Siemens)⁵⁴. Moreover, GABA+ showed greater reliability for both creatine-referenced and water-suppressed estimates^54,55. MM-supressed GABA and GABA+ estimates are also correlated, albeit weakly- to moderately-so^53–55. Consequently, we report GABA+ (henceforth referred to as GABA) to allow for easier comparison of our results to other studies as well as reproducibility.

Basis sets were created using MRSCloud^56,57. For sLASER, ‘localization’ was set to ‘sLASER’, ‘vendor’ to ‘GE’, ‘editing’ to ‘UnEdited’, and ‘TE’ to 35. For MEGA-PRESS, ‘localization’ was set to ‘PRESS’, ‘editing’ to ‘MEGA’, ‘TE’ to 68, ‘edit on’ to ‘1.9’, ‘edit off’ to ‘7.5’, and ‘pulse duration’ to ‘14’. Metabolites included for both basis sets were: Asc, Asp, Cr, CrCH2, EA, GABA, GPC, GSH, Gln, Glu, Gly, H2O, Lac, NAA, NAAG, PCh, PCr, PE, Ser, Tau, mI, and sI. Excitatory-inhibitory ratio (E/I) was calculated as [Glx in i.u.]/[GABA+ in i.u.], a common practice to report E/I using MRS⁵⁸.

Participants were excluded from further analysis if any of their MRS scans had FWHM > 10. Originally we intended to calculate, report and analyze MRS as tissue-corrected quantitative values. However, the Osprey software was reporting Glx values that were an order of magnitude larger than reported Glx values in the literature. Despite several attempts, we were unable to locate the source of this error. Therefore we decided to use Glx/Cr and GABA+/Cr values to be safe.

2.7 Hurst Exponent Calculation

Hurst exponent was calculated from the power spectrum density (PSD) of the BOLD signal. A log-log plot was used, where log power was plotted against log frequency; generally, if a log-log plot results in a linear relationship, it is assumed that the mean slope of this line represents the power-law exponent¹⁸. A PSD shows the distribution of signal variance (‘power’) across frequencies. Complex signals are classified into two categories: fractional Gaussian noise (fGn) and fractional Brownian motion (fBm)^59,60. The former is a stationary signal (i.e., does not vary over time), while the latter is non-stationary with stationary increments⁶⁰. Most physiological signals consist of fBm, but fMRI BOLD is typically conceptualized as fGn once motion-corrected; otherwise put, unprocessed BOLD signal is initially fBm which is converted to fGn with appropriate processing⁶¹. fBm and fGn require distinct H calculation methods⁶⁰. PSD was estimated using Welch’s method⁶² from the Python Scipy.Signal library⁶³. Data were divided into 8 windows of 50% overlap and averaged. The spectral index, , was calculated from the full frequency spectrum. The spectral index was then converted to H using the following equation^60,64:

\[ H = \frac{1 + \beta}{2} \]

Since it cannot be assumed that all fBm is removed from the signal, we use the concept of ‘extended Hurst’ (H’) in this study: for 0 < H < 1, the signal is understood as fGn, while for 1 < H < 2, the signal is assumed to be fBm⁶⁵. More generally, it is also assumed that when 0.5 < H < 1.5, the signal displays 1/f behaviour¹⁸.H was calculated for all voxels in the brain of each subject. A brain mask was then applied which included only GM and the region of the MRS voxel in the visual cortex. H was averaged across the brain mask area, using only non-zero voxels.

2.8 Statistics

All statistical analyses were performed using R⁶⁶ and RStudio (v2023.06.0+421). Glutamate, GABA, E/I ratio, and H were compared between rest and movie conditions using paired Student’s t-tests⁶⁷. Glutamate, GABA, and E/I ratio were tested for correlation with H using Pearson’s method⁶⁸.

3 Results

3.1 Participant Demographics

Twenty-seven participants were originally recruited for the study. Twenty-six of these participants were successfully scanned, but one participant experienced anxiety and chose not to continue. Of the remaining 26 participants, 19 were included in the final analysis: two were removed due to low MRS quality (FWHM > 10) and five were removed due to low fMRI quality (mean FD > 0.15 mm). See Figure 2.

**Figure** 2. **Flowchart of participant, recruitment, scanning and exclusion**. The number of participants at each stage is indicated within each node. Participants were excluded based on MRS FWHM and fMRI Mean FD thresholds, resulting in 19 participants included in the final analysis

The final study sample included 9 males and 10 females between ages 21.3 and 53.4, with a mean age and standard deviation of 30.1 ± 8.7 years.

3.2 Data Quality

3.2.1 MRS

After exclusion, FWHM (mean ± sd) at rest in sLASER and MEGAPRESS were 8.52 ± 0.7 and 7.02 ± 0.7, respectively. During movie watching, FWHM in sLASER and MEGAPRESS were 8.37 ± 0.65 and 6.98 ± 0.86.

Glx and GABA+ were tested for associations with FWHM values during rest and movie watching. Glx during rest was found to be negatively correlated with FWHM (r = -0.53, p = 0.02). No other correlations were significant.

An average of all MRS voxel placements can be seen in Figure 3 A, and a sample of the Osprey sLASER and MEGAPRESS spectrum fits at rest can be seen in Figure 3 B and C, respectively.

In [6]:

layout <- "
  AAAA
  BBCC
  DDDD
  DDDD
  "

mrsvoxel <- image_read("./images/MRSVoxelAverage.png") # From Figures.ipynb
mrsvoxel <- ggplot() + background_image(mrsvoxel) + theme_void()

ospreyslaser <- image_read("../Data_from_BCCHRI/derivatives/slaser/sub-Pilot10/rest/Reports/reportFigures/sub-Pilot10/sub-Pilot10_metab_A_model.jpg")
ospreyslaser_dim <- image_info(ospreyslaser)
ospreyslaser <- ggplot() + background_image(ospreyslaser) + theme_void() #+
  # theme(aspect.ratio = ospreyslaser_dim$height / ospreyslaser_dim$width)

ospreymegapress <- image_read("../Data_from_BCCHRI/derivatives/megapress/sub-Pilot10/rest_osprey/Reports/reportFigures/sub-Pilot10/sub-Pilot10_metab_diff1_model.jpg")
ospreymegapress_dim <- image_info(ospreymegapress)
ospreymegapress <- ggplot() + background_image(ospreymegapress) + theme_void() #+
  # theme(aspect.ratio = ospreymegapress_dim$height / ospreymegapress_dim$width)

# p1 <- free(mrsvoxel) + ospreyslaser + ospreymegapress + glu_FWHM_plot
p1 <- free(mrsvoxel) + free(ospreyslaser) + free(ospreymegapress) + free(gluandgaba_FWHM_plot) +
  plot_layout(design = layout)
p1 + plot_annotation(tag_levels = "A") & theme(plot.tag = element_text(size = 11))

`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'

**Figure** 3. **MRS data quality** A) Average MRS location. Overlay in red of the average voxel location across all nineteen participants, from 0 (white) to 1 (dark red). Darker red represents voxel locations shared by all participants, while white represents voxel locations unique to participants. MRS voxels were registered to MNI space and averaged. Underlay is T1w MNI152 at 0.5 mm resolution. B) Sample Osprey sLASER Spectrum. Yellow spectrum indicates overall model fit. Model fits for individual metabolites are shown in blue below overall fit. C) Sample Osprey MEGAPRESS Spectrum. Yellow spectrum indicates overall model fit. Model fits for individual metabolites are shown in blue below overall fit. D) Correlation plots of Glx / tCr vs. FWHM for rest (green) and movie (red). E) Correlation plots of GABA+ / tCr vs. FWHM for rest (green) and movie (red).

3.2.2 fMRI and Hurst

A sample of the combined grey-matter and MRS voxel mask used to average H values can be seen in Figure 4 A, along with a sample Hurst exponent map in Figure 4 B, and sample fits for H calculation during rest in Figure 4 C.

Mean FD was not correlated with H during rest (r = -0.33, p = 0.16 but was moderately negatively correlated with H during movie watching (r = -0.50, p = 0.03; see Figure 4 D).

In [7]:

# layout="
#   AAAAAAAAA#
#   BBBBBBBBBB
#   CCCCDDDDDD
#   "
layout = "
  AAA
  BBB
  CDD
  "

fmrimask <- image_read("./images/SampleMask.png")
fmrimask_dim <- image_info(fmrimask)
fmrimask <- ggplot() + background_image(fmrimask) + theme_void() #+ theme(aspect.ratio = fmrimask_dim$height / fmrimask_dim$width)

hurstsamp <- image_read("./images/HurstSample.png")
hurstsamp_dim <- image_info(hurstsamp)
hurstsamp <- ggplot() + background_image(hurstsamp) + theme_void() #+ 
  #theme(aspect.ratio = hurstsamp_dim$height / hurstsamp_dim$width)
  #coord_fixed()

psdplot <- image_read("./images/PSDs.png")
psdplot_dim <- image_info(psdplot)
psdplot <- ggplot() + background_image(psdplot) + theme_void() #+ theme(aspect.ratio = psdplot_dim$height / psdplot_dim$width)

meanfd_hurst <- image_read("./notebooks/MeanFD_Hurst_Plot.png")
meanfd_hurst_dim <- image_info(meanfd_hurst)
meanfd_hurst <- ggplot() + background_image(meanfd_hurst) + theme_void() #+ theme(aspect.ratio = meanfd_hurst_dim$height / meanfd_hurst_dim$width)

fmrimask + hurstsamp + psdplot + meanfd_hurst +
  plot_layout(design = layout) +
  plot_annotation(tag_levels = 'A') & theme(plot.tag = element_text(size = 11))

**Figure** 4. **fMRI data quality.** A) sample fMRI grey-matter mask within the MRS voxel. Background: a sample coronal, sagittal, and axial slice is displayed of the mean fMRI scan from the rest acquisition. Foreground: the greymatter/MRS mask used to calculate mean H. B) Sample H map for whole brain. Coronal, sagittal, and axial views of are shown, colour-coded by H values. Colour-coding shows evident tissue differentiation by H value (i.e., gray matter (red), white matter (light red), and cerebrospinal fluid (white)). C) Sample PSD spectrum. All participants’ PSD spectrums during rest are plotted using light grey lines. Mean PSD is plotted in solid red Mean linear regression line is plotted in dashed blue. H for each participant was calculated from the slope of their mean linear regression line. D) Correlation plots of H vs. meanFD for rest (green) and movie (orange).

3.3 Paired Student’s T Tests

Mean ± sd of metabolites, E/I, and H during rest and movie are reported in Table 3. Neither Glx nor GAB+ were different between movie and rest conditions (Figure 5 A & B, respectively). E/I ratio did not change between conditions either (Figure 5 C). H was found to be greater during movie watching than rest (Figure 5 D).

Table 3. Summary of main results

	Rest	Movie	p-value
Glx / tCr	1.19 ± 0.11	1.14 ± 0.11	0.09
GABA+ / tCr	0.45 ± 0.06	0.44 ± 0.06	0.50
E/I Ratio	2.68 ± 0.49	2.65 ± 0.45	0.82
H	0.98 ± 0.98	1.05 ± 1.05	< 0.01

In [8]:

layout <- "
  AABB
  CCDD
  "

glu_boxplot + gaba_boxplot + ei_boxplot + hurst_boxplot + plot_layout(design = layout) + plot_annotation(tag_levels = 'A') & theme(plot.tag = element_text(size = 11))

**Figure** 5. **Paired comparison of metabolite values during rest (green) and movie (orange) conditions.** A) Glx / tCr; B) GABA+ / tCr; and C) E/I. Paired dots represent the same participant across conditions. Mean difference with 95% confidence intervals, and as Cohen’s D are reported at the top of each plot.

3.4 Correlations

H was not found to correlate with Glx, GABA+, or E/I, during rest or movie (Figure 6).

In [9]:

layout <- "
AABB
#CC#
"
glu_hurst_plot + gaba_hurst_plot + ei_hurst_plot + plot_layout(design=layout, guides = 'collect', axes = 'collect', axis_titles = 'collect') + plot_annotation(tag_levels = 'A') & theme(plot.tag = element_text(size = 11))

`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'

**Figure** 6. **Scatter plots of H vs. metabolites.** A) Glx / tCr; B) GABA+ / tCr; and C) E/I. Rest is in green, while movie is in orange. p-values are reported at the top of each plot.

4 Discussion

We report here the first in vivo human study of the E/I-H relationship. An increase in H was observed during movie-watching as compared to rest, indicating a decrease in BOLD signalling complexity in response to visual stimuli. However, no difference in E/I was observed between conditions, nor was H found to be related to E/I during either movie-watching or rest. Our finding that H increases during movie-watching compared to rest is consistent with a previous study by our lab, which found that H increases in the visual network during movie-watching⁶⁵. However, previous studies have observed a task-dependent H decrease using highly-structured, active tasks requiring patient input^69–72. The results of our study in combination with our lab’s previous study⁶⁵, suggest that the naturalistic, passive task of movie-watching induces a different effect on H than more structured tasks. This is in line with literature suggesting that distinct neural responses and BOLD signal characteristics are observed in conventional, active visual tasks compared to naturalistic and passive visual stimuli^65,73.

It may also be that richer scaling properties (higher H) during movie-watching comes about to support the continuous perception of visual stimuli during movie-watching⁶⁵. Moreover, task-induced H changes have been observed to be dependent on task complexity and novelty⁷⁰. Given this finding, it seems reasonable to suggest that differences in task complexity and novelty may account for discrepancies in task-induced H changes.

We found no task-induced change in GABA+. This finding is consistent with a large meta-analytic review which pooled the results of 49 MRS studies⁵¹. They too reported no task-dependent GABA change in the visual cortex. This may due to the technical difficulties of capturing GABA levels using MRS at 3T, owing to its low concentration in the brain as well as signal overlap with higher abundance metabolites⁵¹. However, several studies have successfully reported changes in GABA at 3T using different paradigms / regions (for example:^74–76; see⁵¹ for more).

We also found no change in Glx between conditions. This is somewhat consistent with the literature – while a meta-analysis of recent Glu/Glx studies has reported a small task-induced increase in Glu/Glx in the visual cortex, this finding was non-specific to visual stimuli, including pain, learning, and motor tasks⁵¹. It is also true that a large number of the studies included were conducted at 7T and thus have greater sensitivity to detect changes in Glu/Glx than the present study. Like GABA, glutamate is difficult to capture due to its similar chemical structure and signal overlaps with other metabolites.

Given that no difference was observed for either glutamate or GABA between conditions, it is unsurprising that E/I does not change either.

In terms of the central question of our research, we did not find a relationship between H and E/I during either movie-watching or rest. This is not entirely surprising given the large disparity of existing findings in this regard, especially with regard to the directionality and linearity of the proposed E/I-Hurst relationship^19–25. The heterogeneity across these E/I-Hurst studies highlights the challenges of studying this phenomenon as well as the complexity of any potential relationship between these two metrics. It may be that, given the heterogeneity of these studies in combination with our data, an E/I-Hurst relationship — should it exist — is perhaps more complex than initially thought, and may depend in part on how the data is collected.

Another reason we may not have found a relationship between E/I and H is that we sampled too narrow a range of profiles in terms of E/I and Hurst. Because our study included only healthy adults, it is possible that the participants did not show enough range in E/I and Hurst values to see a relationship. It is also possible that while an E/I-Hurst relationship exists, it is not observed within the visual cortex. This theory seems plausible given that MRS studies of disrupted E/I, mostly conducted within the context of adult autism spectrum disorder, have found changes in E/I within other brain regions such as the anterior cingulate cortex (ACC), frontal lobe, or temporal lobe²⁹. Moreover, findings with reference to changes in excitatory or inhibitory neurotransmitters within the visual cortex tend to be difficult to capture, perhaps indicating that E/I shows less changes in this region⁵¹.

4.1 Limitations

There are two key limitations to our study which should be acknowledged. First is the small sample size. As this was a first-attempt pilot study, only 26 participants were initially scanned. Once individuals were excluded for poor MRI quality, only 19 participants’ data were analyzed. With such a small sample size, it is difficult to make conclusions about a concept as complex as the E/I-Hurst relationship as well as its relations to criticality. The other key limitation of our study is the low sensitivity of our MRS method — as a result, it is hard to conclusively state whether we found no E/I-Hurst link because it truly does not exist or if we found no E/I-Hurst link because our methods lacked the sensitivity to detect it. A possible solution would be to conduct a future study using ultra-high-field 7T MRI as opposed to conventional 3T MRI.

5 Conclusion

In conclusion, our findings do not support a relationship between H and E/I in the visual cortex either during rest or during movie-watching at 3T in humans. In addition, while we found a task-related change in H, we did not find any changes in glutamate, GABA, or E/I between movie and rest. Comparing our findings to the broader literature, E/I balance may be too subtle to be detected with conventional MRS methods. With regards to the broader E/I-Hurst relationship, we similarly suggest that either this relationship is insufficiently captured with our methods, or that the relationship between these two variables may be more complex than originally envisaged—perhaps they are not directly related, but rather connected through other mediating variables in a non-linear fashion. To our knowledge, this is the first in vivo human study to test for this relationship. It is our hope that as the literature grows, more authors will examine this in vivo relationship with respect to other brain regions and using other methods, and will use the lessons learned in this study to inform their own. Hopefully then it will be possible to corroborate findings to probe the complex relationships that may exist with regards to H and E/I in the human brain.

1.

Ogawa, S. & Lee, T.-M. Magnetic resonance imaging of blood vessels at high fields: In vivo and in vitro measurements and image simulation. Magnetic Resonance in Medicine 16, 9–18 (1990).

2.

Ogawa, S., Lee, T. M., Kay, A. R. & Tank, D. W. Brain magnetic resonance imaging with contrast dependent on blood oxygenation. Proceedings of the National Academy of Sciences 87, 9868–9872 (1990).

3.

Stephan, K. E. & Roebroeck, A. A short history of causal modeling of fMRI data. NeuroImage 62, 856–863 (2012).

4.

Deco, G., Jirsa, V. K. & McIntosh, A. R. Resting brains never rest: Computational insights into potential cognitive architectures. Trends in Neurosciences 36, 268–274 (2013).

5.

Campbell, O. L. & Weber, A. M. Monofractal analysis of functional magnetic resonance imaging: An introductory review. Human Brain Mapping 43, 2693–2706 (2022).

6.

Beggs, J. & Timme, N. Being Critical of Criticality in the Brain. Frontiers in Physiology 3, (2012).

7.

Dong, J. et al. Hurst Exponent Analysis of Resting-State fMRI Signal Complexity across the Adult Lifespan. Frontiers in Neuroscience 12, (2018).

8.

Wink, A. M., Bernard, F., Salvador, R., Bullmore, E. & Suckling, J. Age and cholinergic effects on hemodynamics and functional coherence of human hippocampus. Neurobiology of Aging 27, 1395–1404 (2006).

9.

Maxim, V. et al. Fractional Gaussian noise, functional MRI and Alzheimer’s disease. NeuroImage 25, 141–158 (2005).

10.

Long, Z. et al. A Brainnetome Atlas Based Mild Cognitive Impairment Identification Using Hurst Exponent. Frontiers in Aging Neuroscience 10, 103 (2018).

11.

Lai, M.-C. et al. A Shift to Randomness of Brain Oscillations in People with Autism. Biological Psychiatry 68, 1092–1099 (2010).

12.

Dona, O., Hall, G. B. & Noseworthy, M. D. Temporal fractal analysis of the rs-BOLD signal identifies brain abnormalities in autism spectrum disorder. PLOS ONE 12, e0190081 (2017).

13.

Wei, M. et al. Identifying major depressive disorder using Hurst exponent of resting-state brain networks. Psychiatry Research: Neuroimaging 214, 306–312 (2013).

14.

Jing, B. et al. Identifying current and remitted major depressive disorder with the Hurst exponent: A comparative study on two automated anatomical labeling atlases. Oncotarget 8, 90452–90464 (2017).

15.

Dona, O., Noseworthy, M. D., DeMatteo, C. & Connolly, J. F. Fractal Analysis of Brain Blood Oxygenation Level Dependent (BOLD) Signals from Children with Mild Traumatic Brain Injury (mTBI). PLOS ONE 12, e0169647 (2017).

16.

Baranger, M. Chaos, Complexity, and Entropy: A physics talk for non-physicists. New England Complex Systems Institute (2000).

17.

Bassett, D. S. & Gazzaniga, M. S. Understanding complexity in the human brain. Trends in cognitive sciences 15, 200–209 (2011).

18.

Zimmern, V. Why Brain Criticality Is Clinically Relevant: A Scoping Review. Frontiers in Neural Circuits 14, (2020).

19.

Liang, J., Yang, Z. & Zhou, C. Excitation–Inhibition Balance, Neural Criticality, and Activities in Neuronal Circuits. The Neuroscientist 10738584231221766 (2024) doi:10.1177/10738584231221766.

20.

Poil, S.-S., Hardstone, R., Mansvelder, H. D. & Linkenkaer-Hansen, K. Critical-State Dynamics of Avalanches and Oscillations Jointly Emerge from Balanced Excitation/Inhibition in Neuronal Networks. Journal of Neuroscience 32, 9817–9823 (2012).

21.

Lombardi, F., Herrmann, H. J. & de Arcangelis, L. Balance of excitation and inhibition determines 1/f power spectrum in neuronal networks. Chaos: An Interdisciplinary Journal of Nonlinear Science 27, 047402 (2017).

22.

Baumgarten, L. & Bornholdt, S. Critical excitation-inhibition balance in dense neural networks. Physical Review E 100, 010301 (2019).

23.

Bruining, H. et al. Measurement of excitation-inhibition ratio in autism spectrum disorder using critical brain dynamics. Scientific Reports 10, 9195 (2020).

24.

Trakoshis, S. et al. Intrinsic excitation-inhibition imbalance affects medial prefrontal cortex differently in autistic men versus women. eLife 9, e55684.

25.

Gao, R., Peterson, E. J. & Voytek, B. Inferring synaptic excitation/inhibition balance from field potentials. NeuroImage 158, 70–78 (2017).

26.

Tian, Y. et al. Theoretical foundations of studying criticality in the brain. Network Neuroscience 6, 1148–1185 (2022).

27.

Rubinov, M., Sporns, O., Thivierge, J.-P. & Breakspear, M. Neurobiologically Realistic Determinants of Self-Organized Criticality in Networks of Spiking Neurons. PLOS Computational Biology 7, e1002038 (2011).

28.

Rubenstein, J. L. R. & Merzenich, M. M. Model of autism: Increased ratio of excitation/inhibition in key neural systems. Genes, brain, and behavior 2, 255–267 (2003).

29.

Ajram, L. A. et al. The contribution of [1H] magnetic resonance spectroscopy to the study of excitation-inhibition in autism. Progress in Neuro-Psychopharmacology and Biological Psychiatry 89, 236–244 (2019).

30.

Stanley, J. A. & Raz, N. Functional Magnetic Resonance Spectroscopy: The “New” MRS for Cognitive Neuroscience and Psychiatry Research. Frontiers in Psychiatry 9, 76 (2018).

31.

Cordey, H. Jungles. Our Planet (2019).

32.

Li, X., Morgan, P. S., Ashburner, J., Smith, J. & Rorden, C. The first step for neuroimaging data analysis: DICOM to NIfTI conversion. Journal of Neuroscience Methods 264, 47–56 (2016).

33.

Gorgolewski, K. J. et al. The brain imaging data structure, a format for organizing and describing outputs of neuroimaging experiments. Scientific Data 3, 160044 (2016).

34.

Boré, A., Guay, S., Bedetti, C., Meisler, S. & GuenTher, N. Dcm2Bids. (2023) doi:10.5281/zenodo.8436509.

35.

Tustison, N. J. et al. N4ITK: Improved N3 Bias Correction. IEEE Transactions on Medical Imaging 29, 1310–1320 (2010).

36.

Avants, B. B., Epstein, C. L., Grossman, M. & Gee, J. C. Symmetric diffeomorphic image registration with cross-correlation: Evaluating automated labeling of elderly and neurodegenerative brain. Medical Image Analysis 12, 26–41 (2008).

37.

Gorgolewski, K. J. et al. Nipype: A flexible, lightweight and extensible neuroimaging data processing framework in Python. 0.12.0-Rc1. (2016) doi:10.5281/zenodo.50186.

38.

Zhang, Y., Brady, M. & Smith, S. Segmentation of brain MR images through a hidden Markov random field model and the expectation-maximization algorithm. IEEE Transactions on Medical Imaging 20, 45–57 (2001).

39.

Smith, S. M. et al. Advances in functional and structural MR image analysis and implementation as FSL. NeuroImage 23 Suppl 1, S208–219 (2004).

40.

Dale, A. M., Fischl, B. & Sereno, M. I. Cortical surface-based analysis. I. Segmentation and surface reconstruction. NeuroImage 9, 179–194 (1999).

41.

Klein, A. et al. Mindboggling morphometry of human brains. PLoS computational biology 13, e1005350 (2017).

42.

Esteban, O. et al. fMRIPrep: A robust preprocessing pipeline for functional MRI. Nature Methods 16, 111–116 (2019).

43.

Jenkinson, M., Bannister, P., Brady, M. & Smith, S. Improved optimization for the robust and accurate linear registration and motion correction of brain images. NeuroImage 17, 825–841 (2002).

44.

Cox, R. W. AFNI: Software for Analysis and Visualization of Functional Magnetic Resonance Neuroimages. Computers and Biomedical Research 29, 162–173 (1996).

45.

Posse, S. et al. Enhancement of BOLD-contrast sensitivity by single-shot multi-echo functional MR imaging. Magnetic Resonance in Medicine 42, 87–97 (1999).

46.

DuPre, E. et al. TE-dependent analysis of multi-echo fMRI with *tedana*. Journal of Open Source Software 6, 3669 (2021).

47.

Oeltzschner, G. et al. Osprey: Open-Source Processing, Reconstruction & Estimation of Magnetic Resonance Spectroscopy Data. Journal of neuroscience methods 343, 108827 (2020).

48.

Friston, K., Ashburner, J., Kiebel, S., Nichols, T. & Penny, W. Statistical Parametric Mapping: The Analysis of Functional Brain Images. (Elsevier Ltd., 2007).

49.

Gasparovic, C. et al. Use of tissue water as a concentration reference for proton spectroscopic imaging. Magnetic Resonance in Medicine 55, 1219–1226 (2006).

50.

Provencher, S. W. Automatic quantitation of localized in vivo 1H spectra with LCModel. NMR in biomedicine 14, 260–264 (2001).

51.

Pasanta, D. et al. Functional MRS studies of GABA and glutamate/Glx – A systematic review and meta-analysis. Neuroscience & Biobehavioral Reviews 144, 104940 (2023).

52.

Ramadan, S., Lin, A. & Stanwell, P. Glutamate and Glutamine: A Review of In Vivo MRS in the Human Brain. NMR in biomedicine 26, 10.1002/nbm.3045 (2013).

53.

Harris, A. D., Puts, N. A. J., Barker, P. B. & Edden, R. A. E. Spectral-Editing Measurements of GABA in the Human Brain with and without Macromolecule Suppression. Magnetic resonance in medicine 74, 1523–1529 (2015).

54.

Mikkelsen, M. et al. Big GABA II: Water-referenced edited MR spectroscopy at 25 research sites. NeuroImage 191, 537–548 (2019).

55.

Mikkelsen, M. et al. Big GABA: Edited MR spectroscopy at 24 research sites. NeuroImage 159, 32–45 (2017).

56.

Hui, S. C. N. et al. MRSCloud: A cloud-based MRS tool for basis set simulation. Magnetic Resonance in Medicine 88, 1994–2004 (2022).

57.

Mori, S. et al. MRICloud: Delivering High-Throughput MRI Neuroinformatics as Cloud-Based Software as a Service. Computing in Science & Engineering 18, 21–35 (2016).

58.

Rideaux, R. No balance between glutamate+glutamine and GABA+ in visual or motor cortices of the human brain: A magnetic resonance spectroscopy study. NeuroImage 237, 118191 (2021).

59.

Duff, E. P. et al. The power of spectral density analysis for mapping endogenous BOLD signal fluctuations. Human Brain Mapping 29, 778–790 (2008).

60.

Eke, A., Herman, P., Kocsis, L. & Kozak, L. R. Fractal characterization of complexity in temporal physiological signals. Physiological Measurement 23, R1–38 (2002).

61.

Bullmore, E. et al. Wavelets and functional magnetic resonance imaging of the human brain. NeuroImage 23, S234–S249 (2004).

62.

Welch, P. The use of fast Fourier transform for the estimation of power spectra: A method based on time averaging over short, modified periodograms. IEEE Transactions on Audio and Electroacoustics 15, 70–73 (1967).

63.

Virtanen, P. et al. SciPy 1.0: Fundamental algorithms for scientific computing in Python. Nature Methods 17, 261–272 (2020).

64.

Schaefer, A., Brach, J. S., Perera, S. & Sejdić, E. A comparative analysis of spectral exponent estimation techniques for 1/f\(\beta\) processes with applications to the analysis of stride interval time series. Journal of neuroscience methods 222, 118–130 (2014).

65.

Campbell, O., Vanderwal, T. & Weber, A. M. Fractal-Based Analysis of fMRI BOLD Signal During Naturalistic Viewing Conditions. Frontiers in Physiology 12, (2022).

66.

Team, R. C. R: A language and environment for statistical computing. (2021).

67.

Student. The Probable Error of a Mean. Biometrika 6, 1–25 (1908).

68.

Freedman, D., Pisani, R. & Purves, R. Statistics. vol. 4 (W. W. Norton & Company, New York, 2007).

69.

He, B. J. Scale-Free Properties of the Functional Magnetic Resonance Imaging Signal during Rest and Task. Journal of Neuroscience 31, 13786–13795 (2011).

70.

Churchill, N. W. et al. The suppression of scale-free fMRI brain dynamics across three different sources of effort: Aging, task novelty and task difficulty. Scientific Reports 6, 30895 (2016).

71.

Ciuciu, P., Abry, P. & He, B. J. Interplay between functional connectivity and scale-free dynamics in intrinsic fMRI networks. NeuroImage 95, 248–263 (2014).

72.

Barnes, A., Bullmore, E. T. & Suckling, J. Endogenous Human Brain Dynamics Recover Slowly Following Cognitive Effort. PLoS ONE 4, e6626 (2009).

73.

Hasson, U., Malach, R. & Heeger, D. J. Reliability of cortical activity during natural stimulation. Trends in Cognitive Sciences 14, 40–48 (2010).

74.

Floyer-Lea, A., Wylezinska, M., Kincses, T. & Matthews, P. M. Rapid Modulation of GABA Concentration in Human Sensorimotor Cortex During Motor Learning. Journal of Neurophysiology 95, 1639–1644 (2006).

75.

Sampaio-Baptista, C. et al. Changes in functional connectivity and GABA levels with long-term motor learning. NeuroImage 106, 15–20 (2015).

76.

Stagg, C. J., Bachtiar, V. & Johansen-Berg, H. The Role of GABA in Human Motor Learning. Current Biology 21, 480–484 (2011).

Article Notebook